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In this study the numerical performances of wide and compact fourth order formulation of the 
steady 2-D incompressible Navier-Stokes equations will be investigated and compared with 
each other. The benchmark driven cavity flow problem will be solved using both wide and 
compact fourth order formulations and the numerical performances of both formulations will 
be presented and also the advantages and disadvantages of both formulations will be 
discussed. 



1. INTRODUCTION 

In Computational Fluid Dynamics (CFD) field of study, most second order accurate 
finite difference approaches use three point discretizations. Depending on the flow 
problem, if a higher spatial accuracy is desired then fourth order accuracy can be 
chosen. In order to achieve fourth order spatial accuracy, standard five point 
discretization can be used. With five point discretization, wide fourth order 
formulation of the Navier-Stokes formulation has disadvantages near the boundaries 
such that when a wide fourth order formulation is used the points adjacent to the 
boundaries have to be treated specially. 

Another way of achieving fourth order spatial accuracy is to use compact fourth order 
formulations. Compact fourth order formulations provide fourth order spatial 
accuracy in a 3x3 stencil. The main advantage of this type of formulation is that it 
could be easily used at the points adjacent to the boundary without a complexity. 

In the literature it is possible to find many studies that uses "compact high order finite 
difference approximations". These studies can be categorized into two depending on 
the approache used for obtaining the compact high order finite difference 
approximations. In one category, the compact high order finite difference 
approximations are basically generalizations of the Fade scheme and the studies of 
Lele, Visbal and Gaitonde and Zhong are examples for this type of "compact high 
order" studies. For this category of the compact high order finite difference 
approximations Zingg have compared the numerical efficiency of the non-compact 
and compact finite difference methods. 

In the second category, basically the Taylor series and the derivatives of the 
governing equations are used in obtaining the compact high order finite difference 
approximations. As examples of this type of "compact high order" studies, Spotz and 
Carey [18] and Li et al. [13], Zhang [23], Dennis and Hudson [5], MacKinnon and 
Johnson [14] and Gupta et al. [11] have demonstrated the efficiency of compact 
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fourth order formulations of the streamfunction and vorticity equations for uniform 
grids, and also Ge and Zhang [9] and Spotz and Carey [19] have applied the fourth 
order compact formulation to nonuniform grids. This study is intended to contribute to 
this type of "compact high order" literature. 

Recently Erturk and Gokcol [7] have presented a new fourth order compact 
formulation. The uniqueness of this formulation is that the presented compact fourth 
order formulation of the Navier-Stokes equations are in the same form with the 
Navier-Stokes equations with additional coefficients. In this formulation if these 
coefficients are chosen as zero the Navier-Stokes equations are obtained. Therefore if 
a numerical code is written for the solution of the introduced compact fourth order 
formulation, the obtained numerical solution is spatially fourth order ((9(Ax:'*)) 
accurate to the Navier-Stokes equations. In this code if the coefficients are chosen as 
zero then the solution of the introduced compact fourth order formulation is spatially 
second order (0(Ax^)) accurate to the Navier-Stokes equations. Therefore with this 
formulation, using a single equation one can either obtain a second accurate or a 
fourth order accurate numerical solution of the Navier-Stokes equations just by setting 
up some coefficients. Moreover, using the introduced formulation one can easily 
change an existing second order accurate numerical code to provide fourth order 
accuracy and to do this the only thing have to be done is to add some coefficients into 
the existing second order accurate numerical code. When these coefficients are added 
into a second order code to obtain fourth accurate solutions, the code will run slower 
because of the extra CPU work of evaluating these inserted coefficients. Erturk [8] 
have shown that this extra CPU work is slightly dependent on the numerical method 
used. Erturk [8] also showed that for the driven cavity flow when ADI method is used 
in order to achieve numerical solutions with the same level of convergence (defined as 
convergence to the same maximum percent change in flow variables), a fourth order 
accurate numerical code requires 1.8 times more CPU time than a second order 
accurate numerical code. 

Using the Taylor series expansion, one can easily prove that both compact fourth 
order formulations and five point wide fourth order formulations indeed provide 
fourth order spatial accurate numerical solutions. In the literature, to the best of 
authors' knowledge, there is no study that compares the solutions of compact fourth 
order formulation and the solutions of wide fourth order formulation in terms of 
accuracy and numerical performance. 

The aim of this study is then to solve the steady 2-D incompressible Navier-Stokes 
equations with fourth order accuracy using both compact fourth order formulation and 
wide fourth order formulation. First in order to compare the spatial accuracy of both 
formulations we will consider an analytical test case. Then in order to compare both 
formulations in terms of numerical performances, as a test case we will consider the 
benchmark driven cavity flow for different Reynolds numbers with different grid 
mesh and we will compare the CPU time and number of iterations needed for 
convergence for both formulations. We will also document the maximum allowable 
time step that could be used in both formulations as a sign of stability characteristics 
of the formulations. Finally we will document the advantages and disadvantages of 
each formulation. 
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2. FOURTH ORDER COMPACT FORMULATION 



We consider the steady 2-D incompressible Navier-Stokes equations in 
streamfunction {y/) and vorticity ( <y ) formulation. In non-dimensional form they are 
given as 



s-H 7r = -CO 

1 d^co ^ 1 d^a d(o dy/ dco 



Re dx^ Re dy^ dy dx dx dy 
where x and y are the Cartesian coordinates and Re is the Reynolds number. 
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If we discretize equations (1) and (2) using three point central differencing, the 
following finite difference equations provide second order ((9(Ax^,Aj^)) accurate 
solutions 

¥„ + ¥yy=-(o (3) 
1 1 
Re 

where subscripts denote second order ((9(Ax^,Aj^)) three point central finite 
difference derivative approximations. 



0) H CO =11/ CO —w 0) 



(4) 



As explained in Erturk and Gokcol [7], with using second order ((9(Ax^,Ay^)) three 
point central differencing, the solution of the following finite difference equations are 
fourth order ((9( Av;'*,Aj:^Aj^, Aj"* )) accurate to the Navier-Stokes equations (1) and 
(2). 
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We note that equations (5) and (6) are in the same form with equations (3) and (4) 
except with additional coefficients A, B , C , D , E and F . In equations (5) and (6) 
if the coefficients A, B , C , D , E and F are chosen as zero, we identically obtain 
equations (3) and (4). Therefore the solution of equations (5) and (6)become second 
order accurate {0(Ax^,Ay^)) to the Navier-Stokes equations when the coefficients are 
chosen as zero. On the other hand, if the coefficients A, B , C , D , E and F in 
equations (5) and (6) are calculated as they are defined in equation (7), then the 
solution of these equations ((5) and (6)) are fourth order accurate 
{0( Ax* , Ax^ Ay^ , Ay"^ )) to the Navier-Stokes equations. As demonstrated in Erturk [8], 
equations (5) and (6) provide easy way of obtaining either second order or fourth 
order spatial accurate solutions of the Navier-Stokes simply just by using the 
appropriate coefficients. 



In the limit when Ajc — > and Ay — > , the finite difference derivative 
approximations in equations (5), (6) and (7) could be written as partial derivatives as 
the following 
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As it is described briefly in Erturk and Gokcol [7], the numerical solutions of 
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equations (8) and (9) are fourth order accurate ((9( Aj:*, Ax^Aj^, Aj'' )) to the Navier- 
Stokes equations (1) and (2), strictly provided that second order central discretizations 
((9( Ajc^, Aj^)) are used and also strictly provided that a uniform grid mesh with Ax 
and Aj is used. 

One thing is important to note here that, the Navier-Stokes equations (1) and (2) are a 
subset of equations (8) and (9) such that, when the coefficients A, B , C , D , E and 
F in equations (8) and (9) are chosen as zero, we obtain the Navier-Stokes equations 
(1) and (2). Since equations (8) and (9) and the Navier-Stokes equations (1) and (2) 
are in the same form, any numerical method suitable for the Navier-Stokes equations 
(1) and (2) can easily be applied to solve equations (8) and (9). 

Moreover, if we have an existing code that solve the NS equations (1) and (2) with 
second order accuracy, we can easily alter the existing code by adding the coefficients 
A, B , C , D , E and F defined in equation (10) into the code, then the existing 
second order accurate code will turn into a fourth order accurate code. Therefore a 
single numerical code for the solution of equations (8) and (9) can provide both 
second order and fourth order spatial accuracy just by setting some coefficients. 



3. NUMERICAL SOLUTION METHODOLOGY 



We will numerically solve both the NS equations (1) and (2) and the introduced 
compact fourth order formulation of the Navier-Stokes equations (8) and (9). Both of 
these equation sets are nonlinear and therefore they need to be solved in an iterative 
manner. In order to have an iterative numerical algorithm we assign pseudo time 
derivatives to these equations, such that as an example for equations (8) and (9) we 
obtain the following 
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We solve these equations (11) and (12) in the pseudo time domain until the solutions 
converge to steady state. 

For the solution method in the pseudo time, we decided to use the Alternating 
Direction Implicit (ADI) method. ADI method is a very widely used numerical 
method and in this method a two dimensional problem is solved in two sweeps while 
solving the equation implicitly in one dimension in each sweep. The reader is referred 
to [15], [1] and [4] for details. 



3. 1 Compact Fourth Order Solution Methodology 

In order to obtain solutions of the fourth order compact formulation, when we apply 
the ADI method to solve equations (8) and (9), first we solve the streamfunction 
equation and in order to do this we first solve the following tri-diagonal system in x - 
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direction 
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then we solve the following tri-diagonal system in y -direction 
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After this we solve the vorticity equation and in order to do this we first solve the 
following tri-diagonal system in x -direction 
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then we solve the following tri-diagonal system in y -direction 
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In equations (13), (14), (15) and (16), the superscripts denote the iteration time level, 
5^ and 5^ denote the second derivative three point central finite difference 

operators, and similarly 5^ and 5^ denote the first derivative three point central finite 

difference operators in x - and y -direction respectively, for example 

2Ajc 
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where i and j are the grid index and denote any differentiable quantity. 

We note that we use the Thomas algorithm for the solution of these tri-diagonal 
systems . Using the above numerical procedure defined in equations ((13), (14), (15) 
and (16)) we solve equations (8) and (9) and obtain spatially fourth order accurate 
solutions of the Navier-Stokes equations. 



3.2 Wide Fourth Order Solution Methodology 

In order to obtain solutions of the fourth order wide formulation, when we apply the 
ADI method to solve the NS equations (1) and (2), similarly first we solve the 
streamfunction equation and in order to do this we first solve the following penta- 
diagonal system in x -direction 
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then we solve the following penta-diagonal system in y -direction 
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Then, similarly, we solve the vorticity equation and in order to do this we first solve 



the following penta-diagonal system in x -direction 
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then we solve the following penta-diagonal system in y -direction 
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We note that, in equations (18), (19), (20) and (21) the second and first derivative 
operators 5^^ , S^.^ , 5^ and 5^ denote five point central finite difference operators in 

X - and y -direction respectively, for example 
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where / and j are the grid index and denote any differentiable quantity. 



We note that, in numerical solution of the wide formulation we use the modified 
Thomas algorithm for the solution of the penta-diagonal systems. Following the above 
numerical procedure defined in equations (18), (19), (20) and (21) we obtain spatially 
fourth order accurate solutions of the Navier-Stokes equations. 



3.3 Boundary Conditions 

The compact formulations require a 3 x 3 stencil and therefore this formulation can be 
easily used at the first set of grid points near the boundaries. However, wide 
formulations can not be applied directly to the first set of grid points near the 
boundaries. We are going to compare the numerical solutions of wide and compact 
formulations, therefore in order to isolate the effect of the boundary conditions on 
both formulations we would like to use the same boundary conditions in both cases, 
so that the effect of the boundary conditions will be the same in both cases. 

Stortkuhl et al. [21] have presented an analytical asymptotic solution near the comers 

of the cavity and using finite element bilinear shape functions they also have 

presented a singularity removed boundary condition for vorticity at the corner grid 

points as well as at the wall grid points. In this study, we follow Stortkuhl et al. [21] 

and calculate the vorticity values at the wall grid points (circle points in Figure 1-a) as 

the following 

9V 3 1 
= - ^ - -^TTi (V^d + ¥, + ¥,)--i^o).+2co^ + co^ + Ao)^ + a},) (23) 
2Ah 2Ah 8 

where V is the speed of the wall which is equal to 1 for the moving top wall and 

equal to for the three stationary walls. The grid points a, b , c , d , e and / are 

shown in Figure 1-b. 
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For comer points, we again follow Stortkuhl et al. [21] and calculate the vorticity 
values at the comer grid points (diamond point in Figure 1-a) as the following 
9V 3 1 

a>u = ^i//f — {2m + ft>, + 2o)) (24) 

' 2Ah Ah^ ' 4 ' ' ' 

where again V is equal to 1 for the upper two corners and it is equal to for the 

bottom two comers. The grid points b , c, e and / are also shown in Figure 1-c. The 

reader is referred to Stortkuhl et al. [21] for details on the boundary conditions 

In this study, for the first set of grid points adjacent to the wall we decided to use the 
Computational Boundary Method (CBM). For details on the CRM, the reader is 
referred to Huang [12], Yang [22], Gresho [10] and Spotz [20]. Using a fourth order 
one sided approximation for the velocity on the wall we obtain the following 
expression 

^ -25^0 + 48^1 - 36y/^ + I6i//,- 3i//, ^ ^ 
\2h 



w ^ 

on 



where denotes the wall velocity, n is the normal wall direction, h is the grid 

spacing, denotes the grid points on the wall, 1 denotes the first set of grid points 
adjacent to the wall and similarly 2, 3 and 4 denotes the second, third and fourth set of 
grid points adjacent to the wall respectively. From this relation we can calculate the 
streamfunction at the grid points near the boundary (rectangle points in Figure 1-a) as 
the following 

^ + 12/?y, + 25^0 + 36^, - 16^3 + 3^4 
'48 

In order to calculate the vorticity values at these grid points, we used the 
streamfunction equation (1). Using a five point wide fourth order formulation for the 
streamfunction equation, the vorticity at the first set of grid points adjacent to the 
boundary (rectangle points in Figure 1-a) is calculated as the following 

= 

^-^.....^16^,,, -30^ +16^,,, ^ . . 

I2A3;' 

We note that, this approximation require the streamfunction values at grid points 
outside the computational domain (hexagon points in Figure 1-a). In order to find the 
streamfunction values at these grid points we use the following fourth order relation 
-?>Y-x - lOn + 18^1 - 6^^ + ^3 ^ ^ ^^4. ^28) 
I2h 



dn 

where -1 denotes the grid point outside the computational domain. From this relation 
we calculate the value of streamfunction at the exterior grid points ( y/_^ ) and use it in 

equation (27) to calculate the vorticity at the grid points adjacent to the boundary. 

At the first diagonal grid points near the comers (triangle point in Figure 1-a), we 
calculate the streamfunction and vorticity values, using Gauss-Seidel method on the 
five point fourth order wide formulation of the streamfunction equation (1) and 
vorticity equations (2) respectively. 
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We note that in both wide and compact formulations, we used the above described 
boundary conditions for the grid points on the wall and for the first set of grid points 
adjacent to the wall. Then we used the wide and compact formulations to obtain the 
numerical solutions at interior grid points shown in Figure 1-a. 



4. RESULTS AND DISCUSSIONS 

Following the numerical procedure described in the previous section we obtain 
spatially fourth order solutions of the Navier-Stokes equations obtained using both 
compact and wide formulations. 

For the choice of time steps in solving the governing equations, we decided to use 
different time steps, Af, for streamfunction and vorticity equations, as it was also 
done in Erturk [8]. If we solve the streamfunction and vorticity equations using three 
point second order accurate discretization using the ADI method, tri-diagonal matrices 
appear on the implicit LHS of the equations. In streamfunction equation the diagonal 
elements on the LHS matrices become 1+^ , and in vorticity equations the diagonal 

elements on the LHS matrices become 1 + -J-^ . We decided to choose different time 

steps for streamfunction and vorticity equations such that these different time steps 

would make the diagonal elements the same. Therefore we use ^t = alSh^ for the 

streamfunction equation and similarly we use - ccRelsh^ for the vorticity equation, 
where is a coefficient we can choose. We solve both the wide and compact 
formulations using the same time steps. 

In this study we would also like to document the maximum allowable time step 
that we can use in both formulation for convergence as a sign of numerical stability 
characteristics of both formulations. In order to this, we first solve the wide and 
compact formulations using the same time steps. Then we increase ^ , i.e. increase 

the time step , and solve the wide and compact formulations again. We continue 

doing this until at some the solution does not converge. Therefore we would 

document the maximum allowable for convergence for wide and compact 
formulation for a given Reynolds number and grid mesh. 

In this study, in all of the cases considered, we start the iterations from a homogenous 
initial guess and continue until a certain condition of convergence is satisfied. As the 
measure of the convergence to the same level, one option is to use the residual of the 
equations as it was also used by Erturk, Corke and Gokcol [6]. However, we are 
actually solving two different equations, the NS equations (1) and (2) and also the 
fourth order equations (8) and (9), and we are trying to compare the numerical 
performance of these two different equations. Therefore the residual of these two 
different equations may not indicate the convergence to the same level. Alternatively, 
we can use the difference of the streamfunction and vorticity variables between two 
time steps as the measure of convergence for the two equations. However, the 
solutions of the two different equations are slightly different. Since the solutions are 
different, the difference of the streamfunction and vorticity variables between two 
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time steps may not show the same convergence level for these equations also. 
Therefore in this study, as convergence criteria we decided to use the difference of the 
streamfunction and vorticity variables between two time steps normalized by the 
previous value of the corresponding variable, such that 



Residual^^ = max 
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(29) 



These residuals provide an indication of the maximum percent change in i// and co 
variables in each iteration step and in this study we let the iterations converge until 
both Residual^ and Residual^ are less than 10"** . 

Using the Taylor series expansion, mathematically it is straight forward to show that 
the first leading truncation error term in both fourth order compact formulations and 
five point fourth order wide formulations is indeed fourth order. In order to document 
this numerically, we decided to compare the numerical results with a known analytical 
solution using a model problem introduced by Richards and Crane [16]. Inside the 
domain (x, y) = ([0,l],[0,l]) the following analytical solutions satisfy the Navier- 
Stokes equations (1) and (2). 



-e 



Re 

CO = 2e<"^^' (30) 
For this model problem, as the boundary conditions we decided to use the analytical 
solutions defined in equation (30) at the grid points on the boundaries and at the first 
set of grid points adjacent to the boundaries. This way we would be able to avoid any 
effect of numerical boundary condition approximations on the numerical solution and 
concentrate on the accuracy of the solution of both formulations in the interior domain 
for a given analytical values at the boundaries. 

We note that, in equation (30) the vorticity is independent of Re and the solution of 
the streamfunction at different Re numbers looks almost the same in a contour plot, 
therefore for this model problem we only considered the case where Re = 1000 . We 
solve this model problem with wide and compact fourth order formulations using 
different grid mesh with Ah = 1/16 , 1/32 , 1/64 , 1/128 where Ax = Ay = Ah. In these 
fourth order solutions the average absolute difference between the exact solution 
(Wex^'^ex^ given in equation (30) and the numerical solution (y/jj, a>. . ), defined as 

TuVexi^i^yj)- 



E... = ■ 



0), 
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E = 
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(31) 



where A'^ being the total number of grid points, should be proportional to h'" where 
the power should be equal to m = 4 . Since these average absolute differences should 
have the following form 

E = CAh"' (32) 
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logarithm of both sides give 

log£ = logC + OTlogA/i (33) 

In Table 1 we have tabulated the average absolute differences {E^ and E^) for wide 

and compact fourth order solutions for different grid mesh and also in Figure 2, these 
average absolute differences are plotted with respect to the grid spacing in a log-log 
scale. In Figure 2 the average absolute differences of the wide and compact 
formulations are shown by square and circle symbols respectively. In order to have an 
idea on how these average absolute differences change with respect to the grid 
spacing, two linear lines with blue and red colors with slopes of 4 are drawn in the 
same figure. From Figure 2 we can clearly see that both wide and compact fourth 
order formulations (shown by square and circle in the figure), indeed provide fourth 
order accurate solutions, such that for both formulations the slope between log E and 
log Ml is very close to m^A. 

Using the average absolute differences {E^ and E^) tabulated in Table 1, since the 

grid size is decreased by a factor of 2, we can calculate the convergence rate m using 
the following formula 
E{Ah = Ahi) _ 

2 . 

The convergence rate m for both wide and fourth order formulations is also tabulated 
in Table 1 . From this table, we can again see that when the grid spacing is decreased 
progressively by half , the convergence rates of both formulations is very close to 
m « 4 . 



(34) 



Next, we considered the benchmark driven cavity flow problem. For the driven cavity 
flow we obtain fourth order numerical solutions using both wide and compact fourth 
order formulations using the boundary conditions described in the previous section 
and also using both 128 x 128 and 256x256 grid mesh. 

Figures 4, 5 and 6 show the streamfunction contours of the compact and wide 
formulation solutions for a variety of Reynolds numbers (Re = 100, 1000 , 2500) 
obtained with using 256x256 grid points. 

From these figures we can see that both compact and wide fourth order formulation 
solutions are very close to each other. The difference between them is most visible in 
the region of the center of the main circulation. In the literature, for the benchmark 
driven cavity flow, the value of the streamfunction and the vorticity at the center of 
the main circulation and also the location of this center is widely accepted as the flow 
parameters to compare the accuracy of the numerical solutions. In Table 2 we 
tabulated the value of the streamfunction and the vorticity at the center of the main 
circulation and the location of the center for the Reynolds numbers considered. In 
Table 2 we also have tabulated solutions found in the literature that we believe to be 
very accurate. Erturk and Gokcol [7] have presented compact fourth order (Ah*) 
solutions obtained using a very fine grid mesh (600x600). Botella and Peyret [3] have 
used a Chebyshev collocation method and obtained highly accurate spectral solutions 
for the driven cavity flow. Barragy and Carey [2] have used a p-type finite element 
scheme and presented highly accurate (A/z** order) solutions. Schreiber and Keller 
[17], have presented high-order (Ah^ order in theory) solutions obtained using 
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repeated Richardson extrapolation using the solutions obtained on different grid mesh 
sizes. Erturk, Corke and Gokcol [6] have also used repeated Richardson extrapolation 
and presented high-order (A/j^ order in theory) solutions. From Table 2 we can see 
that the compact and wide formulation solutions are very close to each other and also 
all of the solutions agree well with the highly accurate solutions found in the 
literature. 

In order to test the numerical performances of the wide and compact formulations, we 
then solve the benchmark driven cavity flow problem several times with different 
time steps by changing a . In order to see the effect of the number of grids on the 
numerical performances of both formulations, we have tabulated the numerical 
performances of wide and compact formulations to achieve same level of 
convergence. Table 3 shows the CPU time and the number of iterations necessary for 
convergence of the wide and compact fourth order formulations for a 128 x 128 grid 
mesh and also Table 4 shows the same for a 256x256 grid mesh. 

From Tables 3 and 4 we can see that, for the same Reynolds number and for the same 
time step (a) both wide fourth order formulation and compact fourth order 
formulation converge to the same level in about the same number of iteration such 
that the ratio of the number of iteration necessary for convergence for wide 
formulation to the number of iteration necessary for convergence for compact 
formulation is approximately « 1 . For wide and compact fourth order formulations, 
the convergence rate in the pseudo time is approximately the same. 

From Tables 3 and 4 we also can see that, for the same Reynolds number and for the 
same time step (a), the CPU time necessary for convergence of compact fourth order 
formulation is higher than the CPU time necessary for convergence of wide fourth 
order formulation. For the solution of the compact formulation, in each time step the 
coefficients A, B , C , D, E and F in equations (5) and (6) have to be calculated as 
they are defined in equation (7) and this would require a CPU time. The compact 
formulation requires the solution of tri-diagonal systems, and these tri-diagonal 
systems are solved efficiently using the Thomas algorithm. On the other hand, the 
wide formulation requires the solution of penta-diagonal systems. For these penta- 
diagonal systems we used the modified Thomas algorithm. The modified Thomas 
algorithm for penta-diagonal systems runs slower than the Thomas algorithm for tri- 
diagonal systems since it requires more mathematical instructions. Even though the 
modified Thomas algorithm runs slower, from Tables 3 and 4 we see that, for 
convergence the compact formulation requires more CPU time than that of the wide 
formulation. This shows that the extra CPU time necessary to calculate the 
coefficients A, B , C , D , E and F is higher than the extra CPU time necessary for 
the modified Thomas algorithm. For the same Reynolds number and for the same time 
step (a ), in terms of the CPU time necessary for convergence, the wide fourth order 
formulation is more advantageous than the compact forth order formulation such that 
the compact fourth order code requires almost « 1 .5 more CPU time than the wide 
fourth order code. 

For a numerical formulation, the largest time step that the numerical code would not 
diverge is an indicative of the numerical stability of the numerical formulation. In 
order to find the largest time step for the formulations, we progressively increased the 
time step (a) with 0.01 and run the code several times for a given Reynolds number 
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until the solution is no longer convergent such that we could not obtain a solution, i.e. 
the solution is divergent. In Tables 3 and 4 we can see that, in all cases, we can obtain 
numerical solution using larger time steps in compact formulation compared to that of 
wide formulation. This would indicate that the compact formulation is numerically 
more stable than the wide formulation and allows us to use larger time steps. 

When running a code, one would like to use the maximum possible time step (At) for 
a faster convergence. In Table 3 and 4, comparing the CPU time of the compact 
formulation when maximum time step is considered (ar = 1.3) with the CPU time of 
the wide formulation when maximum time step is considered (or = 1.1), we can see 
that when maximum possible time steps are used, the compact formulation runs 
approximately «1.3 times slower than the wide formulation in terms of the CPU 
time. 

2. CONCLUSIONS 

In this study we have numerically solved the Navier-Stokes equations using both 
fourth order compact formulation and fourth order wide formulation and compared 
the numerical performances of both fourth order formulations. 

Solving a model problem which has an exact analytical solution to N-S equations with 
several different grid mesh showed that the solution of both wide and compact 
formulations, indeed converge with fourth degree with respect to the grid spacing. 
Also solving the benchmark driven cavity flow showed that the numerical solution of 
both formulations are very close to each other and produce spatially very accurate 
results. 

We see that both wide and compact formulations converge to a specified convergence 
level in almost the same number of iterations for a given Reynolds number and time 
step. In terms of number of iterations necessary for convergence, both formulations 
have the similar convergence rate. 

For a given Reynolds number and time step, wide formulation requires less CPU time 
than the compact formulation. In terms of computing time, wide formulation runs 
faster and is advantageous over the compact formulation. 

In compact formulation we were able to obtain convergence using larger time steps 
than the time steps we use in wide formulation. This would show that the compact 

formulation is numerically more stable than the wide formulation, i.e. allows us to use 
larger time steps. In terms of numerical stability compact formulation is advantageous 
over wide formulation. 

We would like to point out that in order to use the wide formulation, the points 
adjacent to boundaries have to be treated specially and doing this adds a great 
complexity in coding stage. On the other hand, even though compact formulation does 
not have this complexity and can be easily used for the points adjacent to boundaries, 
in order to use a compact formulation one has to deal with the extra coefficients in the 
equations and this can be counted as the complexity of the compact formulation in 
coding stage. Therefore we can fairly say that both wide and compact formulation 
requires almost the same level of coding effort. 
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Figure 1) Grid points at boundaries 




Figure 2) Average absolute errors in streamfunction and vorticity variables in fourth 
order wide and compact formulations with respect to grid spacing 




Figure 3) Fourth order streamfunction contours of driven cavity flow, i?e = 100 
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Figure 4) Fourth order streamf unction contours of driven cavity flow, Re = 1000 
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Figure 5) Fourth order streamfunction contours of driven cavity flow, Re = 2500 





No. of 

points 




m 




m 


compact 
fomiulation 


16x16 


6.8849E-09 




1 .3766E-08 




32x32 


4.5836E-10 


3.91 


9.1590E-10 


3.91 


64x64 


2.9545E-1 1 


3.96 


5.9001 E-11 


3.96 


128x128 


1.9523E-12 


3.92 


3.8990E-12 


3.92 


256x256 


1.2396E-13 


3.98 


2.4502E-13 


3.99 


wide 
fomiulation 


16x16 


2.641 OE-08 




5.0351 E-08 




32x32 


1 .7978E-09 


3.88 


3.5445E-09 


3.83 


64x64 


1.1710E-10 


3.94 


2.3164E-10 


3.94 


128x128 


7.7605E-12 


3.92 


1.5395E-11 


3.91 


256x256 


4.8750E-13 


3.99 


9.781 3E-1 3 


3.98 



Table 1) Average absolute error and convergence rate of fourth order formulations 



Re 


Ref. 


¥ 


(O 




y 




Present - Compact solution 


-0.1212883 


-1.973684 


0.5195 


0.5430 




Present - Wide solution 


-0.1212327 


-1 .972586 


0.5195 


0.5430 


2500 


Present - Richardson extrap. sol. 


-0.1214392 


-1 .975666 


0.5195 


0.5430 


Erturk and Gokcol [7] 


-0.121472 


-1.976132 


0.5200 


0.5433 




Barragy and Carey [2] 


-0.1214621 




0.5189 


0.5434 




Erturk, Corke and Gokcol [6] 


-0.121470 


-1.976117 








Present - Compact solution 


-0.1188756 


-2.066955 


0.5313 


0.5664 




Present - Wide solution 


-0.1188513 


-2.066586 


0.5313 


0.5664 




Present - Richardson extrap. sol. 


-0.1189358 


-2.067704 


0.5313 


0.5664 


1000 


Erturk and Gokcol [7] 


-0.118938 


-2.067760 


0.5300 


0.5650 


Botella and Peyret [3] 


-0.1189366 


-2.067753 


0.5308 


0.5652 




Schreiber and Keller [1 6] 


-0.118821 


-2.0677 








Erturk, Corke and Gokcol [6] 


-0.118942 


-2.067213 








Barragy and Carey [2] 


-0.11893 










Present - Compact solution 


-0.1035173 


-3.181031 


0.6172 


0.7383 


100 


Present - Wide solution 


-0.1035161 


-3.181007 


0.6172 


0.7383 




Present - Richardson extrap. sol. 


-0.1035190 


-3.181046 


0.6172 


0.7383 



Table 2) Properties at the center of the main circulating eddy in driven cavity flow 







Compact formulation 


Wide formulation 








Re 


a 


compact 


Iteration no,„„p^^, 




Iteration no^^^ 


CPU 

compact 


Ttpr^^tion no 

compact 


CPU 

compact 






Iteration no^^,^ 


CPU,,,, 


max Ar 




1.3 


341 .2 


1 5937 














1.2 


or* o H 

368.1 


1 721 1 














1.1 


398.9 


H 0"7H C 

18715 


OCA O 

256.3 


H O y| OO 

1 8482 


1.56 


1.01 






1.0 


A on n 


OACH "7 

2051 / 


ooo A 

283.0 


OAOcn 
20269 


1.55 


1.01 




2500 


0.9 


4ob. / 


00"7H O 

22/1 2 


O H >l A 

314.0 


oovi a A 
22464 


1.55 


1.01 


1.33 






0.8 


546.6 


25446 


O C H A 

351 .9 


OCO A A 

25200 


1.55 


1.01 






0.7 


625.1 


28944 


O AO O 

398.8 


OOTA A 

28706 


1.57 


1.01 






0.6 


719.4 


33577 


A ^O yl 

463.4 


OOO c o 

33358 


1.55 


1.01 






0.5 


ooo.U 


4000/ 


CCA O 

004.8 


3993/ 


1.55 


1.00 






1.3 


326.7 


1 6529 














1.2 


OCA 

350. 6 


H "7H nc 

17195 














1.1 


oon ^ 


1o/U4 


occ o 


H AA^ O 

1 9018 


1.46 


0.98 






1.0 


A 0"7 C 


2050/ 


OAH n 

291 .9 


OAQCO 

20853 


1.46 


0.98 




1000 


0.9 


yl "7/1 O 

4/4. o 




ooo H 

322.1 


OO AOC 

230ob 


1.47 


0.98 


1.23 






0.8 


C 0"7 A 

527.4 


25430 


OAA O 

360.2 


OC O A H 

25861 


1.46 


0.98 






0.7 


598.3 


28915 


/I H A O 

410.3 


OA A A"7 

29407 


1.46 


0.98 






0.6 


694.0 


oocoo 

33523 


y|"70 O 

478.3 


O A A AC 

34095 


1.45 


0.98 






0.5 


ooo r\ 

833.0 


>l AA >l O 

40048 


568.1 


A ACO"7 

40587 


1.47 


0.99 






1.3 


242.4 


1 1246 














1.2 


OC C "7 

255.7 


H OAOO 

12083 














1.1 


0"7C A 

275.4 


1 31 72 


185.5 


H OACO 

13052 


1.48 


1.01 






1.0 




•\ AA~7Q 
I 44 /O 






1.49 


1.01 




100 


0.9 


336.2 


16071 


224.8 


15883 


1.50 


1.01 


1.31 






0.8 


379.0 


18059 


251.7 


17832 


1.51 


1.01 






0.7 


431.8 


20610 


287.5 


20329 


1.50 


1.01 






0.6 


497.5 


24003 


336.1 


23646 


1.48 


1.02 






0.5 


596.5 


28739 


401.8 


28264 


1.48 


1.02 





Table 3) Numerical performances of compact and wide fourth order formulations with 128 x 128 grid mesh 







Compact formulation 


Wide formulation 








Re 


a 


compact 


Iteration no,„„p^^, 




Iteration no^^^ 


CPU 

compact 


Ttpr^^tion no 

compact 


CPU 

compact 






Iteration no^^,^ 


CPU,,,, 


max Ar 




1.3 


61 56.6 


65039 














1.2 


r^~7 r\r\ c 

6700.5 


69776 














1.1 


7253.1 


76038 


A Cf^C O 

4525.2 


"70H OA 

72189 


1.60 


1.05 






1.0 


/ybi .0 


OOOO/ 


4y/o.o 


/yiys 


1.60 


1.05 




2500 


0.9 


0"700 C 


y2DO/ 


CCAA O 

boOO.o 


o//2o 


1.59 


1.06 


1.36 






0.8 


9770.7 


1 041 02 


6076.6 


AOOO A 

98334 


1.61 


1.06 






0.7 


10998.9 


H H 0"70"7 

1 1 8737 


6834.0 


H H -1 OAO 

111 898 


1.61 


1.06 






0.6 


12838.1 


138195 


OACO A 

8053.0 


H OAOCC 

1 29855 


1.59 


1.06 






0.5 


1 bbUo.l 


1 b5o£:5 


CilO A O 

y/84.8 


1 54981 


1.60 


1.07 






1.3 


5894.5 


62087 














1.2 


6210.2 


67136 














1.1 


5054.0 


54731 


OCOO o 

3523.8 


CCC AO 

55593 


1.43 


0.98 






1.0 


5454. o 


boi oy 


o"7n^ A 

o/yi .0 


CH ACA 

bl4o4 


1.44 


0.98 




1000 


0.9 


bObo.U 


bbbbo 


41 yb.1 


b/bo9 


1.44 


0.99 


1.29 






0.8 


6781 .8 


74835 


A A"7 >1 O 

4674.3 


"7 A AO H 

76031 


1.45 


0.98 






0.7 


"7"7AO "7 

7703.7 


O CO AO 

85302 


CO H A 

5316.6 


86771 


1.45 


0.98 






0.6 


8941 .5 


99195 


^OH H C 

621 1 .5 


HAH r\cc 

101055 


1.44 


0.98 






0.5 


10/47.7 


H H OCOO 

1 1 8538 


"7 >i O yl C 

7434.5 


H OAAOO 

1 20988 


1.45 


0.98 






1.3 


4534.2 


46785 














1.2 


4758.7 


C AOO A 

50239 














1.1 


5054.0 


C y1 "70 ^ 

54731 


o coo o 

3523.8 


C C C AO 

55593 


1.43 


0.98 






1.0 


0404.C5 


OUT uy 


^^/y^ .u 


Dl 404 


1.44 


0.98 




100 


0.9 


6063.0 


66665 


4196.1 


67659 


1.44 


0.99 


1.28 






0.8 


6801.8 


74835 


4674.3 


76031 


1.46 


0.98 






0.7 


7743.7 


85302 


5316.6 


86771 


1.46 


0.98 






0.6 


8941 .5 


99195 


6211.5 


101055 


1.44 


0.98 






0.5 


10747.7 


118538 


7434.5 


120988 


1.45 


0.98 





Table 4) Numerical performances of compact and wide fourth order formulations with 256x256 grid mesh 



